Conservation genomics of urban populations of Streamside Salamander (Ambystoma barbouri)

In Tennessee, populations of the state endangered Streamside Salamander (Ambystoma barbouri) are in decline as their distribution lies mostly within rapidly developing areas in the Nashville Basin. Information regarding the partitioning of genetic variation among populations of A. barbouri and the taxonomic status of these populations relative to northern populations and their congener, the Small-mouthed Salamander (A. texanum), have important implications for management and conservation of this species. Here we combined mitochondrial sequencing and genome-wide single nucleotide polymorphism (SNP) data generated using Genotyping-by-Sequencing (GBS) to investigate patterns of genetic variation within Tennessee populations of A. barbouri, to assess their relationship to populations in Kentucky, and to examine their phylogenetic relationship to the closely related A. texanum. Results from phylogenetic reconstructions reveal a complex history of Tennessee A. barbouri populations with regards to northern populations, unisexual A. barbouri, and A. texanum. Patterns of mitochondrial sequence variation suggest that A. barbouri may have originated within Tennessee and expanded north multiple times into Kentucky, Ohio, Indiana, and West Virginia. Phylogenetic reconstructions based on genome-wide SNP data contradict results based on mitochondrial DNA and correspond to geographic and taxonomic boundaries. Variation in allele frequencies at SNP genotypes, as identified by multivariate analyses and Bayesian assignment tests, identified three evolutionary significant units (ESUs) for A. barbouri within Tennessee. Collectively, these results emphasize the need for prioritizing conservation needs for Tennessee populations of A. barbouri to ensure the long-term persistence of this species.


Introduction
Genetic variation, population structure, and demographic history are increasingly recognized as important factors for the design of effective conservation strategies [1]. Streamside a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Salamander (Ambystoma barbouri, Fig 1) populations in Middle Tennessee are declining due to rapid urbanization in and around the Nashville Basin and as a result, were reclassified from "deemed in need of management" to state endangered in 2018 by the Tennessee Wildlife Resources Agency (TWRA) [2][3][4]. Population fragmentation and loss of genetic variation that inevitably accompany loss of habitat threaten the long-term adaptive potential and persistence of A. barbouri, but can be mitigated by management actions aimed at maintaining genetic diversity. Information regarding A. barbouri's taxonomy, population structure and patterns of genetic variation are currently needed to prioritize conservation needs and to efficiently allocate management resources. Genomic tools are increasingly being used to improve recovery and management planning in at-risk species, including salamanders, by informing taxonomic relationships, demographic histories, and biologically meaningful units of conservation that will preserve genetic diversity and the long-term adaptive potential of species [5][6][7].
Taxonomic uncertainties surrounding the evolutionary relationship between A. barbouri and its sister-species A. texanum (Small-mouthed Salamander) are a concern for precisely identifying targets of conservation. Adults of these two species are nearly indistinguishable based on external morphological features alone and were previously considered to be conspecific. They can be differentiated using scanning electron microscopy on the number and shape of maxillary and premaxillary teeth, larval life-history traits, and choice of breeding environment [8,9]. In addition to their morphological ambiguity, phylogenetic reconstructions using mitochondrial and nuclear sequence data also produced conflicting results regarding the relationship between A. texanum and A. barbouri [10]. Reconstructions based on mitochondrial sequences show that A. texanum and A. barbouri are not reciprocally monophyletic; A. texanum is recovered as a clade nested within A. barbouri from Tennessee, suggesting that pond-breeding A. texanum were more recently derived from a stream-breeding A. barbouri ancestor from Central Tennessee. However, reconstructions based on two nuclear gene sequences resulted in reciprocally monophyletic clades for the two species, consistent with a history where A. barbouri and A. texanum are older independent lineages derived from a shared common ancestor. Additional sequence data are needed to resolve unanswered phylogenetic questions as these earlier reconstructions were based on a limited number of molecular markers (two nuclear and one mitochondrial gene) and representative outgroups.
Mitochondrial evidence also informs on the origins of unisexual ambystomatid populations (family Ambystomatidae) that are common in the Great Lakes region of North America. Unisexual ambystomatids exhibit a unique reproductive mode whereby male sperm activates egg development, but contributes variable amounts of nuclear genetic material depending on compatibility with unrelated cytoplasmic DNA (termed kleptogenesis) [11]. While the nuclear genome of unisexuals is a mixture of different ambystomatid species, all known mitochondrial haplotypes across their range nest within A. barbouri. Mitochondrial haplotypes are most closely related to mitochondrial haplotypes found south of the Ohio River, and southwest of the Kentucky River [12]; however, the relationship of Tennessee A. barbouri to unisexual ambystomatids has not been explored.
At a finer scale, uncertainties also exist regarding the relationship among isolated populations of A. barbouri in Tennessee and discontinuous populations in Kentucky, Indiana, Ohio, and West Virginia. Phylogenetic reconstructions of A. barbouri populations based on mitochondrial sequence data (913 bp) from the D-loop and an adjacent intergenic spacer suggest that Tennessee A. barbouri populations are both ecologically and genetically distinct from more northern populations [13]. At hatching, A. barbouri from Tennessee were smaller and less developed than individuals from Kentucky. Also, laboratory behavioral assays show that Tennessee A. barbouri were similar to western Kentucky individuals, but less active than individuals from northern populations [10]. These differences were supported by genetic data; phylogenetic reconstructions recovered mitochondrial haplotypes from Tennessee as monophyletic and basal to haplotypes from Kentucky, Ohio, and West Virginia [13]. These results support an early divergence of Tennessee populations and raise questions regarding the geographic origin of A. barbouri. However, conclusions from this study were limited as these reconstructions are based on a single mitochondrial gene from only three individuals from a few populations within a kilometer of each other in Rutherford County.
Observations from field surveys suggest that A. barbouri populations in Tennessee are in decline, and the accompanying loss of genetic variation in small populations further threatens the long-term persistence of this species [14]. Maintenance of genetic variation is a fundamental priority for conservation planning and requires information regarding the partitioning of genetic variation within and between isolated populations. Data on patterns of genetic variation in Tennessee populations of A. barbouri are limited to a handful of mitochondrial and nuclear sequences used for phylogenetic studies [13,10]. Genome-wide surveys of genetic variation from fine-scale population sampling across the state are critical for assessing genetic variation within populations, establishing units of conservation and maintaining historical patterns of gene flow between populations to improve the outcome of recovery efforts [15].
Here we investigate patterns of genetic variation within and among populations of Tennessee A. barbouri as well as the taxonomic relationship between A. barbouri and A. texanum in order to prioritize conservation needs and inform management practices aimed at maximizing long-term persistence of this species. Mitochondrial sequence data and genome wide SNP genotypes were used to (1) further examine. the phylogenetic relationship between A. barbouri and A. texanum, (2) review the taxonomic relationship of disjunct populations of A. barbouri in Tennessee, relative to northern populations, (3) evaluate patterns of genetic differentiation between geographically isolated populations of A. barbouri within the state of Tennessee and (4) estimate within-population genetic variation and examine demographic history of Tennessee A. barbouri. Taxonomic relationships, units of conservation, and geographic partitioning of genetic variation are discussed in the context of establishing conservation priorities and designing effective management strategies for this species.

Tissue collection and DNA extraction
Tissue samples of A. barbouri in the form of adult tail clips, eggs, and whole larvae were obtained from field surveys or from collaborators (Fig 2, Table 1). Whole larvae were euthanized with a solution of tricaine sulphonate (buffered MS-222) at a concentration of 3g/L. This protocol was approved by the Institutional Animal Care and Use Committee at Tennessee Tech University (protocol 17-18-006P) and scientific collection permits were obtained from Tennessee Wildlife Resources Agency (permit 1475). Historical and predicted sites of A. barbouri were surveyed from January 2018 to March 2018 and again from November 2018 to March 2019, coinciding with timing of oviposition as reported by Niemiller et al. [8]. At each site, A. barbouri adults and eggs were collected by turning over cover objects within and near seasonal streams. Additionally, pools were searched for free-swimming larvae; only one larva was sampled from any single pool to avoid sampling related individuals from the same clutch. A total of 235 individuals were included for genetic analysis. Samples included 225 individuals from 13 populations of A. barbouri in Tennessee spanning six counties in the Nashville Basin: Bedford County (B6), Davidson County (D3), Rutherford County (R1, R7, and R9), Sumner County (S2, S5, S7, and S8), Wilson County (W1, W3, and W4), and Williamson County (Wil2). A total of five A. texanum individuals were sampled; these included four individuals from (Arnold Airforce Base, Coffee County, TN, AAFB) and one from Craighead Co., Arkansas (AKT1). Outgroups included two A. mabeei (North Carolina), two A. maculatum (AAFB), and one A. talpoideum (AAFB).
DNA was extracted from tail clips, eggs, and larvae using the EZNA Tissue DNA Mini kit (OMEGA BIO-TEK) following the manufacturer's protocol except that DNA was eluted in water. Approximately 1200 bp of the mitochondrial D-loop was targeted for PCR amplification using primers developed by Shaffer and McKnight [15] ( Table 2). The mitochondrial D-loop has been shown to be informative for evaluating population structure and species relationships within the genus Ambystoma [12,[15][16][17][18]. Polymerase chain reactions (PCR) included 2 μl extracted DNA, 0.5 μM of each primer, and 1X GoTaq Master Mix (Promega) to a total volume of 20 μl. Conditions for PCR were as follows: initial denaturation step of 2 minutes at

Phylogenetic reconstructions for mitochondrial haplotypes
Phylogenetic reconstructions were estimated for all unique mitochondrial D-loop haplotypes and an additional 53 A. barbouri, 83 unisexual ambystomatids, and 27 A. texanum sequences obtained from Genbank. Both maximum likelihood (ML) and Bayesian optimality criteria were used for phylogenetic analyses. Maximum-likelihood analyses were performed using the software RAxML [21] on the CIPRES Science Gateway [22] under the GTR+G model. Nodal support was estimated using rapid bootstrapping (1000 replicates). Bayesian phylogenetic reconstructions were performed using MrBayes 3.2.1 [23] also on the CIPRES Science Gateway. The best model of substitution was selected by Modeltest [24], implemented in MEGA X [25] using Bayesian information criterion (BIC). The Markov chain Monte Carlo (MCMC) algorithm ran for 10,000,000 generations, sampling every 1,000 generations. Two independent runs were performed and the resulting trees were combined after the deletion of a burnin (25%). Convergence of the MCMC chain was assessed by estimating the effective sample size for paramaters (>200) and by visualizing parameter trace plots using the software Tracer 1.7.1 [26,27]. A majority-rule consensus tree was generated and nodal support was estimated by posterior probabilities.

GBS library preparation and sequencing
A total of 227 individuals were included for GBS sequencing including 201 individuals from 12 populations of A. barbouri in Tennessee, 17 individuals from five populations of A. barbouri in Kentucky, five individuals from two populations of A. texanum, and two sampled outgroups A. talpoideum (N = 1) and A. mabeei (N = 1; Table 1). Genotyping-by-Sequencing libraries were prepared using the restriction enzyme ApeKI following the protocol of Elshire et al. [28]. Genomic DNA was quantified using Quant-iT Picogreen dsDNA Assay Kit (Thermo Fisher Scientific), and all samples were standardized to between 5-6.5 ng/uL (50-65 ng total genomic DNA). Extracted DNA was digested with the restriction enzyme ApeKI. Adaptors containing PCR binding sites and individual barcodes were ligated onto digested DNA. Barcoded DNA was pooled and PCR amplified using primers that bind to the ligated adaptors (see Elshire et al. [28] for primer sequences). The resulting PCR products were cleaned using the Qiagen PCR purification kit and then cleaned again using the AxyPrep Mag PCR Clean-up kit Table 2. Primers used for PCR amplification of the mitochondrial D-loop. Primers THR (forward) and 651 (reverse) were used for initial amplification of entire~1300 bp. Internal primers 007 (forward), DL3 (reverse), and DL1 (reverse) were used for Sanger sequencing (Shaffer & McKnight, 1996).

SNP discovery and filtering
The Stacks program process_radtags was used to filter and demultiplex raw reads based on barcoded sequences [29]. We used the seven-step de novo clustering pipeline ipyrad v. 3.5 [30] to generate and filter SNP datasets used in downstream analyses. Quality filtering of raw sequence reads converted bases with Phred scores <33 to Ns, reads with more than 5 Ns were removed. Reads were clustered using a sequence similarity threshold of 90% both within and between sampled individuals, with a minimum read depth of six. Individuals with fewer than 500,000 reads were excluded from downstream analyses. Loci with observed heterozygosity (H o ) greater than 0.5 were removed to filter out possible paralogs. The final SNP dataset was then filtered to remove loci deviating from Hardy-Weinberg equilibrium (p < 0.05), loci genotyped in less than 60% of individuals, and SNPs with a minor allele frequency (MAF) less than 0.01. Only one SNP was retained per sequenced tag.

Population genetic variation and effective population size
Estimates of within population genetic variation for 12 Tennessee A. barbouri populations were estimated using the R-package DiveRsity v. 1.9.9 [31]. Summary statistics were estimated separately for each population and included the proportion of SNPs that were polymorphic within each population (P), allelic richness (A R ), and observed and expected heterozygosity (H o and H e ). Estimates of P, A R , H o , and H e were estimated for the entire SNP dataset and for the subset of SNPs that were polymorphic within each population. Effective population size (N E ) was estimated for Tennessee A. barbouri populations using both the linkage disequilibrium (LD) method and the heterozygote-excess method [32] as implemented in the program N E Estimator v2 [33]. Populations with sample sizes of less than 20 individuals were excluded from analyses as these datasets did not provide enough signal for reliable N E estimation [34]. The minor allele frequency (MAF) parameter was set at 0.01 and 95% confidence intervals were estimated using the parametric chi-squared method.
Population structure. Population pairwise F ST values [35] and G" ST [36] were estimated for the 12 Tennessee A. barbouri populations using the R package diveRsity [31]. Hierarchical partitioning of genetic variation across Tennessee populations of A. barbouri was examined using Analysis of Molecular Variances (AMOVA) as implemented in Arlequin 3.5 [37,38]. Results from phylogenetic reconstructions based on mitochondrial D-loop haplotypes were used to generate hypotheses regarding higher-level structuring of populations. Significance of variance components was determined using 1000 permutations.
The optimal number of genetic clusters (K) based on genomic SNPs was estimated using both a multivariate approach and a Bayesian-based assignment method. Discriminant Analysis of Principal Component (DAPC) was performed using the 'adegenet' package in R (v. 2.1.4) [39,40]. The optimal number of principal components (PCs) was determined by cross-validation using the 'xvalDapc' function (1000 replicates) and the number of PCs with the highest assignment success was retained for DAPC analyses. The find.clusters function was first used to identify the optimal value of K based on a Bayesian Information Criterion (BIC) process. The optimal K was then used to perform a DAPC analysis to describe the relationship between the genetic clusters. Individual membership probabilities were assessed using the function compoplot. DAPC analysis was performed on three hierarchical datasets; these included 1) 12 populations of A. barbouri in Tennessee, 2) all sampled Tennessee and Kentucky A. barbouri populations, and 3) all A. barbouri populations and four A. texanum individuals from AAFB.
Bayesian assignment tests were performed on Tennessee A. barbouri populations using the program STRUCTURE 2.3.4 [41]. Kentucky A. barbouri populations were not included in assignment analyses due to small sample sizes at these sites. Values of K (number of populations) ranged from 1 to 12 populations with 20 replicate runs per value of K; MCMC simulations were performed for a burn-in of 250,000 iterations and an additional 1,000,000 iterations were retained for the final analysis. Results were summarized using the software package CLUMPAK [42]. The optimal number of groups (K) was determined using the Delta-K method [43] as implemented in STRUCTURE HARVESTER v0.6.94 application [44]. The Distruct (1.1) application in CLUMPAK was used to produce the final barplots.

Multispecies coalescent using BPP
A multispecies coalescent (MSC) model was used to examine taxonomic boundaries, estimate species trees, and estimate population divergence times (t = years before present) as implemented in the software BPP v. 4.4 [45]. Input datasets were formatted as full-length sequence alignments from our GBS sequencing and included all A. barbouri and A. texanum populations. Loci represented by fewer than 100 individuals were removed from the analysis. Joint species-delimitation and species-tree reconstructions were performed using the A11 model, where individuals were assigned to sampled populations. Run parameters were as follows: uniform rooted trees were used as the species model prior, the theta prior (θ = 4N E μ) was assigned a gamma distribution with α = 3 and β = 0.04 and the prior for tau (τ = tμ) was gamma distributed with α = 3 and β = 0.2. Three independent MCMC runs were implemented, each with a burn-in of 5,000, a sample frequency of 10 and a total of 50,000 iterations. Species delimitation and species-tree topologies with the highest posterior probabilities from the A11 runs were used as input for the A00 analyses. The A00 model estimates divergence time parameters and long-term N E , assuming fixed species assignments and a fixed species tree. The priors for A00 were the same as priors used for the A11 model. Three independent runs were also performed for the A00 analysis, each with a burn-in of 5,000, a sample frequency of 10 and a total of 50,000 iterations. MCMC chains were pooled and absolute estimates of τ were converted to years before present (YBP) in the BPPR using the genome wide mutation rate estimate for the vertebrate nuclear genome where μ = 1.21 x10− 9 [46,47].

Mitochondrial sequencing
Sanger sequencing of the mitochondrial D-loop resulted in an~1100 bp sequence alignment in the 81 individuals selected for sequencing. A total of 60 unique haplotypes were identified from the newly sequenced individuals, 39 haplotypes were identified from A. barbouri populations sampled in Tennessee (See S1 Table for GenBank accessions), and 15 haplotypes were identified from populations sampled in Kentucky. An additional six haplotypes were sequenced from outgroups. The final alignment contained 126 variable sites (excluding outgroups); of these, 120 sites were parsimony informative. No haplotypes were shared across multiple populations.

Phylogenetic reconstructions
Both Bayesian and ML phylogenetic reconstructions of mitochondrial D-loop haplotypes identified six major clades across A. barbouri and A. texanum haplotypes (Fig 3). The Clade I

GBS sequencing
The average number of retained sequence reads per individual generated from sequencing of GBS libraries was 5.86 x 10 6 ± 5.9 x 10 6 ; 220 individuals were retained after removal of low coverage individuals (defined as < 500,000 total reads). A total of 1.56 x 10 6 loci were recovered from the de novo assembly in ipyrad and 440,629 loci were retained after filtering. Additional filtering for SNPs (i.e. filtering for low representation, HWE, MAF, and max heterozygosity) recovered 1,169 SNPs for the dataset that included only Tennessee A. barbouri, 516 SNPs in the dataset that included all sampled A. barbouri populations in Tennessee and Kentucky, and 500 SNPs in the dataset that included all A. barbouri and A. texanum populations.

Within Population Genetic Variation
Much of the genetic variation within our SNP dataset was in the form of fixed differences between populations ( Table 3). The proportion of SNPs that were polymorphic within populations ranged from 0.045 (W3 and Wil2) to 0.267 (S7). When only polymorphic loci were considered, the Will2 population had the highest estimates of H o and H e ; however, this population had a low proportion of polymorphic SNPs overall. Rutherford County populations (R1 and R7) had the lowest estimates of H o and H e , when considering only polymorphic loci and for the entire SNP dataset.

Effective population size
Estimates of effective population size (N E ) were generated for five populations with adequate sampling (N�20), including S7, S8, R1, R7, and D3 (Table 4). Results from LD analyses ranged from N E = 15 (R1) to N E = 108 (108) (S8); the upper confidence interval for S8 was 1, indicating low signal for this dataset. Results from the heterozygote excess method were inconclusive (1) for two of the five analyzed populations. Estimates from the remaining three populations ranged from 24 (D3) to 406 (S7); however, CIs included 1 for two of these estimates (S7 and S8).

Population structure
Pairwise F ST estimates for A. barbouri populations within Tennessee were geographically structured in a hierarchical fashion (Table 5) The three monophyletic clades of Tennessee A. barbouri recovered from mitochondrial phylogenetic reconstructions were further examined using genome-wide SNP genotypes in an AMOVA framework (Table 6). Most variance in SNP genotypes could be attributed to differences between the three clades (63.64%). Differences between populations within clades only Table 3. Standard measures of genetic diversity for 12 populations of Ambystoma barbouri in Tennessee based 586 SNP loci. Summary statistics include the average number of individuals genotyped per locus (N) and the proportion of SNPs that were polymorphic within each population (P). Allelic richness (A R ), observed heterozygosity (H O ), and expected heterozygosity (H e ) were calculated for all loci and again for only those loci that were polymorphic within each population. Results of BIC for DAPC analysis suggested that three genetic clusters (K = 3) best represented genetic variation among populations sampled within Tennessee (Fig 4A, S1 Fig). Membership of each cluster was geographically partitioned and consisted of a northern cluster that included populations from Sumner County, a central cluster comprised of individuals from Wilson, and Davidson counties, and a southern cluster that represented Bedford, Rutherford, and Williamson County populations. Membership probabilities and DAPC plots indicated some degree of admixture between the central and southern clusters, but no admixture within the northern cluster (S1 Fig). A separate analysis that included Tennessee A. barbouri and the five sampled populations from Kentucky also resulted in three genetic clusters (Fig 4B). Kentucky populations of A. barbouri were isolated from populations in Tennessee. Central and southern populations were assigned to the same cluster and populations in Sumner county were assigned to the third cluster. When A. texanum individuals were included in a third analysis, results of the BIC analysis indicated the number of genetic clusters was five (Fig 4C). Tennessee populations were assigned to northern, central, and southern clusters and scatter plots showed Tennessee populations to be tightly grouped with the fourth cluster that included all Kentucky populations. All A. texanum individuals formed the fifth cluster that was well separated from all A. barbouri individuals.

All Loci Polymorphic Loci
Analysis of Bayesian assignment tests for Tennessee populations of A. barbouri were consistent with results from DAPC analysis that suggested the presence of three geographically partitioned genetic clusters within the state (Fig 5). The optimal value of K, as determined by the Delta K method, indicated four genetic clusters; however, increasing K from 2 to 3 did not

Multispecies coalescent analyses
A total of 988 loci were retained for the MSC analyses performed by BPP. All three A11 runs recovered the same taxonomic groups and species tree topology with the highest posterior probability. Out of 18 sampled populations, 14 populations were genetically distinct. Within Tennessee, A. texanum and all A. barbouri populations, with the exception of R1-R7, were identified as distinct groups (Fig 6). Kentucky A. barbouri populations were split into two groups; populations FC, JC, RR, and SL were grouped together and the SW population was its own genetic group. The separation of individual populations in this analysis does not indicate that each population warrants recognition at the species-level as simulations have shown that BPP will split groups at the population-level when many loci are analyzed [48]. The A11 species tree reconstruction with the highest posterior identified the same Tennessee clades recovered in our mitochondrial gene trees with the exception of the placement of A. texanum. In the species tree, A. texanum was basal to all populations of A. barbouri in Tennessee and Kentucky. The estimated divergence time (A00 analysis) between A. texanum and A. barbouri was 1.6 million YBP. All Sumner A. barbouri formed a monophyletic group with Kentucky populations as a monophyletic sister group. Divergence time estimates indicated that Sumner populations and Kentucky populations shared a common ancestor < 500,000 YBP. All A. barbouri populations south of the Cumberland River were monophyletic (TMRCA 260,000 YBP) and were further split into two groups that corresponded to clades

Discussion
Patterns of genomic variation and taxonomic relationships identified here have important implications for developing management strategies aimed at the long-term conservation of A. barbouri. Despite the documented decline of A. barbouri populations in the Nashville Basin, results from GBS-derived SNP genotyping indicate that estimates of genetic variation in extant populations of Tennessee A. barbouri are similar to SNP-based datasets examined from other ambystomatids [34,49]. Partitioning of genetic variation between Tennessee populations suggests that hydrogeography of the Nashville Basin has shaped patterns of gene flow. Both mtDNA and genomic SNP genotypes showed similar patterns with respect to population structure within the state of Tennessee that should be used to inform the designation of conservation units. Our results also reveal a complex history of Tennessee A. barbouri with more northern populations and with its sister species A. texanum. Despite marked differences in phylogenetic reconstructions based on mtDNA sequencing and nuclear genotypes, collectively these results indicate that Tennessee is genetically unique from northern A. barbouri and A. texanum.

Genetic variation within populations
Loss of genetic variation resulting from rapid population decline threatens the long-term success of conservation efforts due to loss of adaptive potential and fixation of deleterious alleles [50]. For A. barbouri populations in the Nashville Basin, maintenance of adaptive genetic variation is particularly critical as these populations are faced with habitat alteration by urbanization in addition to climate change [51]. Evaluating the genetic health of populations is valuable for management and conservation planning; however, the interpretation of genetic variation estimates from SNP data is challenging due to the limited number of comparable studies utilizing reduced-representation methods in salamanders. Heterozygosity estimates obtained here are similar to estimates obtained in the handful of published SNP-based genetic studies. In a recent survey of SNPs from ddRAD sequencing, estimates of H O ranged from 0.165 to 0.269 for populations of the Mole Salamander (A. talpoideum) and were slightly higher for populations of the Marbled Salamander (A. opacum; 0.207-0.298) [37]. Another SNP-based survey examined the effects of land use on genetic variation in the Northern Two-lined Salamander (Eurycea bislineata). Fusco et al. [52] reported nearly identical estimates of H O from urban, suburban and rural salamander populations (H O = 0.265, 0.278 and 0.275, respectively) and concluded that genetic variation had been maintained, despite habitat disturbance. However, estimates of heterozygosity may not be directly comparable across independent studies as filtering parameters in bioinformatic pipelines can influence population genetic summary statistics. Stringent filtering criteria may preferentially retain loci in conserved regions of the genome and downwardly bias estimates of genetic variation [53]. There have been previous population genetic studies in A. texanum and A. barbouri populations based on microsatellite markers. Micheletti and Storfer [54] estimated genetic variation at 11 microsatellite loci in 76 populations of A. barbouri distributed throughout Kentucky, Ohio, and Indiana; the average H e from their study ranged from 0.67-0.81, respectively. Few conclusions can be made by comparing results from these prior microsatellites studies to our results based on SNP genotypes; microsatellites are known to have a much higher rate of mutation than nucleotide substitutions, increasing the expected amount of allelic diversity for a given population size [55].

Effective population size
A reduction in estimates of contemporary effective population size (N E ) may also signal a decline in the genetic health of at-risk populations. Our estimates of N E were low (N E from the LD method averaged 58.0 across five populations), but were not unusual for LD-derived estimates of N E in salamanders. Published estimates of N E reported for amphibians have typically been under 100 [49,56,57]. Effective population size estimates reported here were similar to values reported for the congeneric endangered California Tiger Salamander (A. californiense; N E values of   [58], Leora's Stream Salamander (A. leorae; N E values of   [59] and for the Long-toed Salamander (A. macrodactylum; N E = 23-2070 [60]. In A. macrodactylum, Savage et al. [61] estimated N E for 47 breeding populations and more than half of these estimates were less than 50. Life history factors may provide some explanation for low N E in this group as salamanders often exhibit high variance in reproductive success and larval survival within populations. Variance in reproductive success will increase relatedness among individuals in a population, inflating linkage disequilibrium across loci and reducing N E . Obtaining robust estimates of N E and other demographic parameters in at-risk species can be challenging due to the sample sizes required for LD detection. Many threatened taxa are inherently rare and prohibitively difficult to sample in large numbers. For robust estimation of N E , it is recommended that sample sizes be greater than 30 for most systems [62]. When the number of sampled individuals is much smaller than the effective size, LD-based N E estimates can be downwardly biased and confidence intervals can be large (or infinite) due to inadequate signal in the dataset [63]. We limited our N E analysis to populations with sample sizes > 20 individuals, which left only five populations for N E estimation. The infinite upper confidence interval for our N E estimate in population Sumner 8 was likely due to an inadequate sample size for this population (N = 21).

Phylogenetic reconstructions and population structure
Mitochondrial gene-tree reconstructions suggest a complex biogeographic history for A. barbouri in Tennessee that involves populations of unisexual ambystomatids, A. texanum, and populations of A. barbouri from the northern core region (Kentucky, Indiana, and Ohio). Our results support findings by Bogart et al. [12] demonstrating a monophyletic relationship of unisexual ambystomatid mtDNA haplotypes nested within A. barbouri. The phylogenetic reconstruction shown here further details this relationship whereby unisexual ambystomatids likely shared a maternal common ancestor with A. barbouri populations in Kentucky and populations from Middle Tennessee. This pattern, coupled with phylogenies based on nuclear genes, points to a history of repeated hybridization events initiated by an A. barbouri maternal ancestor in the southern portion of its range.
With regards to their mitochondrial lineage, A. texanum is nested within present-day A. barbouri, rendering A. barbouri paraphyletic. This phylogenetic pattern contradicts the scenario proposed by Kraus and Petranka [9], that stream-dwelling A. barbouri descended from the more widespread pond-breeding A. texanum. Species-tree reconstructions based on genomewide SNP data are incongruent with relationships recovered from the mitochondrial trees. This species-tree topology adheres to conventional taxonomic and geographical boundaries (Fig 6), where A. texanum is a sister-species to all A. barbouri populations. The TMRCA for A. barbouri and A. texanum based on SNP genotypes was estimated at~1.6 million YBP making this a relatively recent split between these ecologically distinct species [64]. Incongruence between mitochondrial gene trees and reconstruction based on nuclear data are not uncommon and are often attributed to ancestral lineage sorting, introgression, and/or sex-biased dispersal [65,66]. Determining the cause of mitonuclear discordance would require additional sampling across the distributions of both A. barbouri and A. texanum; however, it is relevant to note that recent studies have demonstrated that hybridization is common across salamander lineages (including the genus Ambystoma) and may facilitate rapid diversification [67].
The hydrogeography of the Nashville Basin, together with cyclic glacial movements during the Pleistocene, may have shaped contemporary patterns of genetic variation in Tennessee populations of A. barbouri. Mitochondrial-based phylogenetic reconstructions and partitioning of genetic variation at SNP loci (including assignment tests, DAPC analyses, AMOVA, and MSC reconstructions) identify the same three genetically distinct groups in Tennessee; these include a northern cluster, a central cluster, and a southern cluster. The Cumberland River appears to have served as a major barrier to gene flow between the northern and the central/ southern clusters. The central and southern clusters both occur south of the Cumberland River and are divided by smaller regional drainage patterns. Individuals in the northern and central clusters occupy the Old Hickory Lake and Lower Stones River watersheds, respectively. The central cluster is bordered by the Cumberland River to the North and by the East Fork of the Stones River to the east. The southern cluster includes populations from three watersheds including the West Fork Stones River Watershed, Upper Duck River Watershed, and Mill Creek Watershed; the southern cluster (with the exception of the Wil2 population) is bordered by the Stones River to the north and by the Duck River to the south.
Mitochondrial-based phylogenetic reconstructions revealed that Tennessee A. barbouri populations are not monophyletic with respect to northern populations, which may reflect a history of repeated range expansions from Tennessee into the northern end of A. barbouri's present day distribution. The estimated timing of divergence between the three Tennessee clusters falls within the mid-late Pleistocene. Glacial advances during the Pleistocene did not extend into Kentucky and Tennessee and this region likely served as a refuge for both terrestrial and aquatic fauna [68]. As glaciers retreated, populations at the northern limits of these refugia would have advanced northward into newly habitable territory. These expansions would have left a genetic signature characterized by reduced genetic variation in the northern end of their distribution as a result of founder effects. Repeated expansions of A. barbouri northward from Tennessee may explain the polyphyletic relationships of contemporary Tennessee A. barbouri with respect to populations in the north. Specifically, mitochondrial A. barbouri populations in northeastern Kentucky, Ohio, and Indiana (Clade V) are nested within Tennessee A. barbouri as are mitochondrial haplotypes from ambystomatid unisexuals and A. texanum. Genetic evidence of repeated northward expansions has been reported for the Eastern Woodrat (Neotoma floridana) [69] and other southeastern fauna (reviewed in Hewitt 2000; [70,71]).

Management implications and conclusions
Results from our genomic survey have specific implications for the design of management and conservation strategies that may improve the long-term persistence of A. barbouri in Tennessee. First, patterns of genetic variation in mitochondrial and nuclear genomes support the assignment of three genetically distinct units for management that warrant the designation of evolutionary significant units (ESUs), where ESUs are defined as groups of populations that show phylogeographic differentiation for mtDNA haplotypes and divergence in allele frequencies in nuclear markers [72,73]. These three units include a Northern ESU encompassing all Sumner County populations (S2, S5, S7, and S8), a Central ESU that includes populations from eastern Davidson and Wilson Counties (D3, W1, W3 and W4), and a Southern ESU that includes populations in Williamson, Rutherford, and Bedford counties (Wil2, B6, R1, R7, and R9). Patterns of differentiation between these three groups of populations suggest a long history of genetic isolation at both mitochondrial and nuclear markers, such that different groups are likely to possess unique combinations of adaptive genetic variation and have likely experienced independent evolutionary trajectories.
Evolutionarily and phylogenetically distinctive populations contribute disproportionately to genetic diversity and should be ranked highly in regards to conservation priority. Geographically peripheral populations are frequently observed to be genetically less-redundant than more central populations and are more likely to possess potentially adaptive genetic variation [74]. For A. barbouri, Tennessee populations are geographically disjunct from the majority of A. barbouri breeding populations [14] and occupy the southernmost edge of the species distribution. Results from this study suggest that Tennessee A. barbouri should be prioritized for conservation planning as these populations are both genetically diverse and evolutionarily distinct from populations in the northern part of their distribution. The distinctiveness of these populations is further evidenced by observed differences in reproductive life-history traits, including mean diameter of early stage ovum size and number of eggs per clutch [8]. Prioritizing peripheral populations with adaptive genetic variation and evolutionary potential is even more critical when considering environmental challenges that accompany climate change. Amphibians in general are very sensitive to climate change as their reproductive life histories are linked to temperature and precipitation [75]. However, there is evidence that populations at the warm-range edge of their distribution are more resilient [76]. Protection of A. barbouri breeding sites in Tennessee may be instrumental to ensuring the long-term viability of this species as a whole.